Hamiltonian dynamics of homopolymer chain models 
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The Hamiltonian dynamics of chains of nonlinearly coupled particles is numerically investigated 
in two and three dimensions. Simple, off-lattice homopolymer models are used to represent the 
interparticle potentials. Time averages of observables numerically computed along dynamical tra- 
jectories are found to reproduce results given by the statistical mechanics of homopolymer models. 
The dynamical treatment, however, indicates a nontrivial transition between regimes of slow and 
fast phase space mixing. Such a transition is inaccessible to a statistical mechanical treatment and 
reflects a bimodality in the relaxation of time averages to corresponding ensemble averages, ft is also 
found that a change in the energy dependence of the largest Lyapunov exponent indicates the 0- 
transition between filamentary and globular polymer configurations, clearly detecting the transition 
even for a finite number of particles. 
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I. INTRODUCTION 

Many important facets of polymer physics can be ap- 
proached by means of equilibrium statistical mechanics 
0] Dynamical aspects are typically treated by 
means of numerical Monte Carlo simulations or Langevin 
dynamics 0. ■ Recent developments (see Ref . @ for 
a review and an extensive bibliography) in the theory 
of nonlinear Hamiltonian systems with many degrees of 
freedom provide an alternative viewpoint from which to 
approach polymer physics, and are expected to convey 
information complementary to that obtainable by tradi- 
tional methods. It is the aim of the present work to use 
these recent developments to investigate the Hamiltonian 
dynamics of simplified homopolymer models. 

The geometrical interpretation of nonlinear Hamilto- 
nian systems with highly dimensional phase space has 
seen considerable progress in the past decade, both in 
terms of theory and in application to physical models 
H, 0, 0. Concepts from Riemannian geometry have 
been used to map the natural evolution of a system at 
constant energy onto geodesic trajectories of a manifold 
in phase space. This leads to a natural definition of 
Hamiltonian chaos, which is quantified by the largest 
Lyapunov exponent, a measure of how quickly trajec- 
tories diverge when initially separated by an arbitrarily 
small distance in phase space. 

The study of chaoticity as a function of energy density 
(energy per degree of freedom) is a rich topic with physi- 
cal applications. It has been found that in many systems 
there exists a value of the energy density at which a sud- 
den change in the energy dependence of the Lyapunov 
exponent occurs [rH ITsf . For energy densities above 
this value, there is a region of strong chaos, allowing 
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a system to efficiently sample phase space, while below 
this value there is a region of weak chaos, in which the 
system cannot readily sample all available phase space. 
Such a transition signifies a strong stochasticity threshold 
(SST), and is associated with a change in the topology of 
accessible phase space |l9j . This threshold has also been 
associated with conventional phase transitions |20j |. and 
thus provides a novel tool by which to characterize and 
study them. We also note that the efficiency of phase 
space sampling is intimately related to a physical relax- 
ation time, thus providing information about the time 
necessary for a system to reach equilibrium. 

The ability of Hamiltonian dynamics to quantify the 
efficiency of phase space sampling and to detect conven- 
tional phase transitions makes it a natural candidate for 
the study of polymer dynamics, a field replete with stud- 
ies using equilibrium statistical mechanics and computer 
simulation methods to characterize a polymer phase di- 
agram. The Hamiltonian dynamics approach may be 
particularly useful to provide complementary informa- 
tion on the phase space of protein-like polymers, that 
present a cooperative transition from a disordered un- 
folded ensemble to a compact folded state. To complete 
the folding transition, a protein-like polymer needs to 
navigate a complicated, high-dimensional phase space, 
that can be populated with metastable intermediate 
states and can exhibit a certain degree of local rough- 
ness [H IH E H Q EE 13- The purpose of this 
paper is to show the feasibility and utility of applying 
Hamiltonian dynamics to the study of simple homopoly- 
mer models as motivation for the future study of protein 
dynamics. 

The Hamiltonian dynamics of off-lattice homopolymer 
models are studied at three levels of approximation. Af- 
ter discussing the general model in Sec. [HI we consider 
the free chain approximation in Sec. lIIII the self-avoiding 
approximation in Sec. IIVI and the Lennard- Jones model 
in Sec. E] We compare the scaling of the mean squared 
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end-to-end distance computed along dynamical trajecto- 
ries with theoretical values co mput ed using equilibrium 
statistical mechanics 0, EH H2 El EH . We also fo- 
cus on the relaxation time and Lyapunov exponent as a 
function of energy density. For the Lennard- Jones model, 
we find that the Lyapunov exponent gives a signature 
of the nonconventional 0-transition between filamentary 
and globular configurations of the homopolymer. This is 
a strong indication that Hamiltonian dynamics may pro- 
vide a useful tool in characterizing the topolo gy o f the 
accessible phase-space of protein-like polymers [26| . 



II. THE HOMOPOLYMER CHAIN MODEL 

We represent, as it is common in simple models, a poly- 
mer as a chain of N+l beads connected by springs, where 
a configuration of the chain is defined by the positions 
{qo,...,qAr} of the beads in e?-dimensional continuous 
space. The Hamiltonian is defined as 
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where the {pi} are the canonically conjugate momenta. 
A simple choice for the interaction potential p5ll2^.l2^ | 

is 
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where ry = |ry| = |q s : — q_j\ gives the interparticle dis- 
tances. The expression / sp r(?"ij)j representing spring-like 
anharmonic interactions between neighboring beads in 
the homopolymer chain, is given by 
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where tq is the equilibrium distance between nearest 
neighbors along the chain. 

We can simplify d translational and d(d — l)/2 rota- 
tional degrees of freedom by enforcing the constraints of 
vanishing total momentum and angular momentum 
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on the initial condition. This way the total energy E is all 
internal energy and we can interpret the dynamics gen- 
erated by the Hamiltonian as a constant temperature 
evolution. 

The Hamiltonian equations of motion, 



dH 
dpi' 



dH 

dm' 



i = 0,...,N (5) 



are numerically integrated using an efficient third order, 
bilateral, symplectic algorithm |2{j. The initial condi- 
tions are given by random configurations of beads com- 
patible with the constraints (QJ. 

We expect convergence of relevant observables to en- 
semble averages after a transient, nonequilibrated regime. 
Once equilibrium is attained, time averages will only fluc- 
tuate about ensemble averages. We are interested in the 
dependence of the relaxation time tr on the energy per 
degree of polymerization e = E/N. 

Recall that the ergodic hypothesis underlying statisti- 
cal mechanics is not sufficient to ensure the convergence 
of time and ensemble averages in finite time, and that 
such convergence requires a stronger condition of a phase 
space mixing dynamics. Riemannian geometry provides 
a useful framework for the study of Hamiltonian systems 
with many degrees of freedom. This framework has been 
detailed elsewhere [3(]], and is briefly reviewed in App.IO 
of this paper. A necessary condition for a mixing dynam- 
ics is that trajectories in phase space are unstable with 
respect to arbitrarily small variations of the initial con- 
ditions. A quantitative measure of the mean instability 
of nearby trajectories is given by the largest Lyapunov 
exponent Ai (see App. IbI for the definition), which also 
gives a measure of the time scale necessary for the loss 
of memory of initial conditions. 

It has been shown that reg imes of both fast and slow 
phase space mixing exist |l7l Il8l [l9j • The crossover be- 
tween the two is characterized by a sudden change in 
relaxation time tr as a function of energy density e, and 
is intimately connected with a crossover between different 
scaling laws for Ai(e). In the language of previous work, 
this crossover defines a transition from weak to strong 
chaos, which governs the efficiency in which phase space 
is sampled. 



III. THE IDEAL POLYMER 

We initially consider the potential energy U of Eq. © 
with the constants r)ij and a set equal to zero: 
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This defines the simplest idealization of a flexible polymer 
chain, and is the Hamiltonian dynamical counterpart of 
the random walk model (RWM) on a periodic lattice. 
Average properties of the simple lattice RWM are easily 
computed and visualized |iJ- L3- LI ■ 

For the ideal homopolymer model in the large N limit, 
we can analytically compute the e-dependence of the 
largest Lyapunov exponent by taking advantage of a 
method developed in [33l The theoretical predic- 

tion of Ai (e) , together with other statistical averages that 
are analytically obtained for the ideal model, are then 
compared with the time averages of the same quantities 
numerically computed along the dynamical trajectories. 
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Moreover, in the dynamical simulations, we also estimate 
the relaxation time tr needed for convergence of an ob- 
servable - in our case, the mean value of the square end- 
to-end distance defined below - to its expected value. 
The dynamical simulations of the system are performed 
for a finite number of monomer units, N. It is thus also 
interesting to compare the analytic results, obtained in 
the limit N — > oo, with the numeric results obtained at 
finite values of N. 



A. The end-to-end distance 

The size of a polymer chain is usually characterized by 
its end-to-end distance R = |qjy — Q.o|- F°r a RWM, the 
mean value of R 2 can be computed within the canonical 
ensemble (see App. IaI for the details): it results to be 
proportional to N 



(R 2 ) = N£ 2 ((3) . 



(7) 



The coefficient £ is a monotonically decreasing function 
of the inverse temperature (3 = l/^ksT), whose value 
tends to ro as j3 — ► oo. 

In the thermodynamic limit, N — > oo while e is con- 
stant, we can obtain the microcanonical average (A)^ of 
any observable function A({qi}) in the parametric form 

m 



(A)M = (A) G ((3) 



(4(e) 
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where (A) G is the canonical average of the observable A 
and Z(J3) is the canonical partition function. For finite 
N we have (/)„(/?) = (/) G (/3) + 0(±). 

We therefore expect that, for sufficiently large N, and 
after a suitable transition time, the time averages of the 
end-to-end distance (R 2 ) of the chain, numerically com- 
puted along the dynamical trajectories of the system, will 
agree with the scaling law (R 2 ) oc N at any fixed value 
of the energy density e. 

In Fig. n we show the comparison between analytical 
and numerical results for two different values |47| of the 
energy density e. The agreement between ensemble and 
time averages is very good, even though the simulated 
values of N are not very large. 



B. The Lyapunov exponent and the relaxation time 
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FIG. 1: The power law R 2 oc N is correctly recovered in dy- 
namical averages with a off-lattice random walk model. Re- 
sults are shown for two different values of the energy density: 
e = 0.1 (full circles) and e — 1000 (open circles). Solid and 
dotted line represent the best fit lines, entirely consistent with 
the theoretically predicted behavior. 



|33l l34j as a function of e is reviewed in App. [C] This 
method is used to compute Ai for the ideal homopolymer 
chain, the results of which are shown in Fig. [21 

Figure shows that numerically computed values of 
Ai for N = 25, N = 60, and N = 120 are in excel- 
lent agreement with analytical results. The existence of 
a crossover energy density e c between two different scal- 
ing laws of the exponent Ai is also clear. Such a bi- 
modality is characteristic of a SST, at which a transition 
between regimes of a qualitatively different chaoticity oc- 
curs [Til Hsj . The transition from above to below e c is 
coupled with a change in the topology of phase space tra- 
jectories that results in slower phase space mixing |20| . 
To further support this, we have numerically computed 
tr(c), the relaxation time necessary for time averages of 
the mean square radius (R 2 ) to achieve theoretically pre- 
dicted values. It is evident from Fig. that there exists a 
rather sudden change in the behavior of tr(c) at e = e c , 
clear sign of the occurrence of a SST. 



We now explore the dependence of the largest Lya- 
punov exponent Ai on the energy density e in order to 
identify a crossover in the scaling behavior, that is the 
signal of a transition between slow and fast phase space 
mixing (i.e., a transition between slow and fast relax- 
ation of time averages to ensemble averages). A method 
to analytically compute the largest Lyapunov exponent 



IV. THE SELF-AVOIDING POLYMER 

The RWM can be straightforwardly enriched by intro- 
ducing repulsive interactions between monomers. The 
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new potential energy is 
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FIG. 2: Largest Lyapunov exponent Ai vs. energy density e 
for the random walk model: comparison between theoretical 
prediction (solid line) and numerical results (circles); dotted 
lines are references to asymptotic power laws e 2 and e 2 / 3 , and 
their intersection defines the strong stochasticity threshold. 
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FIG. 3: Relaxation time tr (defined as the time of the first 
occurrence of the theoretical value of (R 2 )) vs. energy density 
e for the random walk model. The picture shows the results 
obtained for a chain with N = 22. 
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giving the counterpart to the self-avoiding walk (SAW) 
on a lattice. In lattice models of polymers, short ranged 
repulsive interactions between monomers are incorpo- 
rated through an excluded volume condition, which re- 
quires that any two monomers maintain at least one lat- 
tice space separation. This model, which is more physical 
than the simple RWM, accurately describes the swollen 
phase of a real polymer chain in a good solvent. 

The mean squared end-to-end distance for this model 
scales as (R 2 ) oc N 2v . The existence of a nontrivial ex- 
ponent v is known by renormalization group calculations 
|3fij . but the value, which is expected to depend on di- 
mensionality, is not analytically known. Flory |37| de- 
vised a simple approximation scheme for computing the 
exponent v in any dimension. Basic assumption is that 
N monomers are uniformly distributed within the vol- 
ume R d of the swollen polymer: all interactions between 
monomers are disregarded. Given a monomer, the prob- 
ability that another one will fall into its excluded volume 
v is therefore Nv/R d . Since there are N monomers, the 
repulsive energy is 



rep 



kT 



R d 



(10) 



The increase in the size of the polymer comes together 
with an increase of the entropy, that for an ideal (i.e. 
without volume exclusion) chain can be estimated as 
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(11) 



The free energy F — E rcp —TS is minimized by the radius 
R = i?F given by 



R d+2 



v£ 2 N 3 , 



which implies that 



2 + d 



(12) 



(13) 



In spite of the rather drastic approximations, this is a 
remarkable assessment, as it gives the correct value of 
v for d = 1 and values within a percent of the most 
accurate numerical results J3(| for d = 2 and d = 3. 
However, Flory's theory doesn't work when the degree of 
polymerization N is too small for the approximation of 
non interacting monomers to hold. To be more specific, 
if the repulsive energy in the ideal configuration R = 
Rq ~ N x l 2 l is less than kT, the polymer doesn't swell, 
so that v — 1/2 as in the RWM case. The repulsive 
interactions are relevant only when the so-called chain 
interaction parameter Q 



(14) 
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FIG. 4: Mean squared end-to-end distance for different values 
of the ratio cr/ro (from bottom to top, 0.1, 0.3, 0.5, 0.7, 0.9) for 
the model defined by Eq. @ . The solid lines follow the power 
law of the 3-dimensional self-avoiding walk N 2l/ , the dashed 
lines represents the behavior linear in N of the random walk 
model. 



• < 0.1 gives (R 2 ) oc N (the value extracted from 
the best fit is 0.99 ± 0.02), as expected from the 
RWM; 

• — > 0.5 gives (R 2 ) oc Ni (numerical values from 
the three linear fits are between 1.19 and 1.21 with 
errors of order 0.02), which is the expected SAW 
scaling; 

• — = 0.3 can't be satisfactorily fitted with just one 
line: a change in the scaling of (R 2 ) occurs at an 
intermediate value of N. 

We assume that the information in Fig. (0} can be syn- 
thesized as follows: 



(R 2 ) 
(R 2 ) 



kiN 
koN 2u 



N < N c 
N > N r 



where 



N r = a — 



(16) 



(17) 



with the constant a independent on a or 7'o. Continuity 
for N = N c dictates the following relation between k\ 
and k2'. 



to 
k 2 



4-d 

Q,2 + d 



(18) 



is large enough, z>l. The above formula implies that 
for d > 4 the effect of repulsion becomes negligible and 
the polymer may be considered ideal. 

In the following, we use the framework of Flory's the- 
ory to interpret the scaling behavior of the off-lattice 
SAW for different values of the ratio a/ro in a unified 
way. 



A. Scaling behavior of (R 2 ) 



We can now express N in terms of the parameter z de- 
fined by Eq. I|14[l. thus obtaining the scaling relation 



where 



m = 



4-d 

lor z < a 2 



a d+2 z {d+2){4-d) f or z > a' 



(19) 



(20) 



The repulsive potential \^r-J m t ne Hamiltonian © 
is sufficiently steep that we can assume two monomers i 
and j cannot get closer than a distance r»j ~ a. The ex- 
cluded volume is therefore v « a d , while the equilibrium 
distance between two nearest neighbors along the chain 
is given by £ ~ rrj. In view of Flory's theory, we expect 
to observe a linear dependence on N of (R 2 ) when 



N<N r 



ro 



(15) 



N 2 



while for N > N c we should observe (R 2 

Figure^presents the results of numerical investigation 
of the dynamical behavior of self-avoiding homopolymers 
containing between 10 and 100 monomers for several val- 
ues of f-. The analysis of data shows that: 



Within the approximations of Flory's theory, we can 
therefore treat in a unified way the transition between 
the random and the self-avoiding walk behavior. As a 
consistency check, Fig. shows how the curves of Fig. 01 
collapse onto a unique curve f(z) once expressed in terms 
of z. 



B. Dynamical behavior 

The behavior of the Lyapunov exponent of the self- 
avoiding chain is invest igat ed in this section. Figure [H| 
compares the numerical [38| results for the SAW (cr/r = 
VlO 11 ) and the analytical results for the RWM (cr = 0). 
It is clear that the analytic curve computed with a = 
provides a good approximation for the numerically com- 
puted Lyapunov exponent as a function of energy density, 
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FIG. 5: Collapse plot for the data of Fig. @. The dotted and 
solid line are the two asymptotic behaviors in Eq. 1201 . 




FIG. 6: Comparison between the largest Lyapunov exponent 
for the random walk case (circles) and the self-avoiding walk 
case (triangles). The solid line is the analytic prediction for 
the random walk case. 



even when afro = VlO 11 . Even though the scaling law 
changes from (R 2 ) ~ N for the RWM to (i? 2 ) - N* 
for the SAW, the dynamics, as characterized by Ai(e), 
maintains features of the random walk. That is, the in- 
troduction of the repulsive energy in the Hamiltonian has 
no apparent effect on the instability of dynamical trajec- 
tories in phase space. 

This fact can be explained as follows: Because the re- 
pulsive potential grows very steeply for rij < 1, the phase 
space for the SAW can be thought of as that for the RWM 
with inaccessible regions. If these forbidden regions of 
phase space are not too large (a/r < 1), we expect that 
sufficiently uniform observables in phase space will not 
be significantly affected. In the strongly chaotic regime, 
the average Ricci curvature and its fluctuations, which 
determine the value of Lyapunov exponent (see Adp. IC)| . 
do not change significantly when measured in different 
regions of the ambient manifold. As such, they are not 
sensitive to the exclusion of some regions of phase space, 
and thus we expect the Lyapunov exponent to be negli- 
gibly affected. We expect that this reasoning may cease 
to be valid when the range of the repulsive potential ap- 
proaches the interatomic distance {a/r > 1). 



V. THE LENNARD-JONES POLYMER 

The results from previous sections have demonstrated 
that the application of Hamiltonian dynamics to the 
study of homopolymers is both feasible and provides 
dynamical insight otherwise inaccessible. We now con- 



sider a homopolymer chain in which attractive forces are 
present in addition to the repulsive forces considered in 
the previous section. Recalling our potential J5J), we set 
Vij = where 77 is a constant, for all i and j. For more 
complicated heteropolymers such as minimalist protein 
models, the values of 77.^ can be selected to effectively 
represent the interactions between different amino-acid 
types (see, e.g., ^(|)- This, however, is beyond the scope 
of our current work and is left for future study. 

The intermonomer potential (J2J) can be rewritten in 
the form 



JV 

U L j = C/rwm + 5m 47 

i=l j<i 



(21) 



where 7 = i] 2 /4 and A = a j ^/fj. 

It is known that a homopolymer chain with this po- 
tential can adopt two distinct phases associated with the 
dominance of either the attractive or repulsive interaction 
energy .3?) ■ With the Lennard- Jones potential above, the 
relative importance of the attractive and repulsive ener- 
gies can be varied by varying the temperature. At high 
temperature, the attractive part of the potential is neg- 
ligible and the chain exists in a swollen phase, closely 
mimicking the SAW model in which only repulsive forces 
are relevant. The power law of the SAW, (i? 2 ) oc N 2v , 
is recovered. At low temperatures, the attractive term 
becomes relevant and the chain collapses into a com- 
pact configuration, resulting in the radius of the system 
scaling as (R 2 ) ~ N 2 / d . The separation between the 
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FIG. 7: The Lennard- Jones polymer (I2H in dimension 2 ex- 
hibits different scaling laws above (full circles, e = 10) and 
below (open circles, e = —2) the 0-point. Solid line is refer- 
ence to (R 2 ) oc N 3/2 , dotted hne to (R 2 ) oc N. 



self- avoiding and collapsed chain regimes is marked by a 
0-temperature which indicates a phase transition. We 
study the properties of the Lennard- Jones homopolymer 
model in dimension d = 2, and thus expect Flory's esti- 
mate of v = 3/4 to hold for the SAW, and v = 1/2 to 
hold for the collapsed chain. 



A. Scaling laws 

Dynamical simulations are performed in the range 
N G [10,100]. We fix the parameters 7 and A so that 
at high temperature the system gives the scaling of Sec- 
tion [^^] for this range of N. In the dynamical simula- 
tions we vary the energy density e, which is a function 
of the temperature according to Eq. (JSJ), and in our mi- 
crocanonical simulations acquires the meaning of mean 
energy per degree of freedom. Figure {7\ shows that two 
different scaling laws are found for values of the energy 
density e = 10 and e = — 2, implying that at these values 
of e the system is in different phases. In fact, the value of 
the exponent jumps between v = 3/4 and v = 1/2 when 
the energy density reaches a value eg . This transition has 
been proven to be a tricritical point [lj. 

The proper way to estimate the value eg at which the 
transition occurs, as well as the tricritical exponent ug, 
would be to perform a finite size scaling analysis along 
the lines illustrated in Ref. ^(j- Unfortunately, we can- 
not follow this route because the fluctuations of the end- 



to-end distance in our off-lattice model are much more 
pronounced compared to the lattice SAW models usually 
employed for this measure. A rough estimate of eg can 
be obtained nevertheless by studying the variation of the 
exponent v with respect to the energy density. With the 
choice of parameters a — 0.1,6 = 10, A = 9,7 = 2, the 
tricritical exponent vg — 4/7 |22| corresponds to an en- 
ergy density between —0.8 and —0.7 (the yellow line in 
Fig. EJ|. 



B. The Lyapunov exponent and dynamical 
characterization of the 0-transition 

It has been recently conjectured that, in the presence 
of a second order phase transition, there exists a close 
relationship between dynamical properties described by 
Lyapun ov exponents and statistical ensemble averages 
[2*5 liH lil I43L H3. In particular, the behavior of the 
Lyapunov exponent is expected to change abruptly as a 
function of energy density e near the transition. This 
provides motivation to seek whether there is a dynami- 
cal signature marking the nonconventional 0-transition 
expected in our homopolymer system. 

To test for a dynamical indication of the O-transition, 
we have numerically computed the largest Lyapunov ex- 
ponent Ai (e) of the system as a function of energy den- 
sity e. A sharp change in the slope of Ai(e) takes place 
for a value of the energy density close to the transition 
value eg estimated above, as can be seen in Fig.0 Even 
though the 0-transition is only properly defined in the 
thermodynamic limit N — *■ 00, numerical calculations of 
the Lyapunov exponent for N = 50 (blue points) and 
N = 100 (red points) already appear to indicate the ex- 
istence of a phase transition. 

We have analytically computed the high temperature 
behavior of Ai (e) in order to validate the numerical cal- 
culations by comparing to the behavior of observables in 
the limit N — > 00. Even though the statistical averages 
are not exactly computable, we use the method described 
in App. El in which the attractive part of the potential 
is negligible and the chain is, in practice, a self-avoiding 
one. We note that the behavior of the Lyapunov expo- 
nent for the self-avoiding polymer in three dimensions is 
equivalent to the random walk. 

A similar result holds for the self-avoiding chain in 
two dimensions, but a more careful estimate is necessary 
to compute the average quantities described in App. [U] 
This is because the integrals diverge when integrating 
over all possible configurations using the canonical mea- 
sure exp(— (3H) with the free chain potential C/rwm @- 
Physically, however, it is very unlikely that two particles 
will get closer than the length scale 7 of the Lennard- 
Jones interaction. Thus, we can approximate the differ- 
ence between the canonical measure for the free chain 
and the canonical measure for the self-avoiding chain by 
neglecting a spherical volume of radius 7 around each 
particle. Such a procedure does not change the result in 
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FIG. 8: Lyapunov exponent as a function of the energy den- 
sity for N = 50 (blue) and N — 100 (red) Lennard- Jones 
chains. The purple points show the behavior of a N — 50 
SAW, while the purple line is the theoretically predicted value 
for an infinite length SAW. The yellow line marks the esti- 
mated value of ee- 



three dimensional space because the integrand functions 
vanish in these regions. In Fig. [5] we show the analytical 
result for the Lyapunov exponent obtained by applying 
this method. The agreement with the corresponding nu- 
merical values is very good at energy density greater than 
the critical value eg, which corresponds to the high tem- 
perature phase. At lower energy density, the attractive 
part of the Lennard-Jones potential ceases to be negligi- 
ble, and agreement with the values for the SAW case is 
lost. 



VI. CONCLUSIONS 

We have shown that the Hamiltonian dynamics de- 
scription of continuous space homopolymer chain models 
reproduces the statistically predicted scaling laws, even 
for homopolymer chains with a small number of monomer 
units. This result is nontrivial since the dynamical model 
allows monomer units to move freely in continuous space, 
with changes in conformation being driven by nontrivial 
differentiable and conservative dynamics. The dynamical 
description of homopolymer systems also conveys qual- 
itatively new information, giving a physical time scale, 
tr, that measures the amount of time necessary for time 



averages to converge to their statistical counterparts for 
a system with arbitrary initial conditions. Thus, if tr 
were to exceed the observational time lapse, experimental 
measurements would be inadequate to predict statistical 
mechanical averages. Moreover, tr has nontrivial en- 
ergy density dependence, and a sudden change of tr(e) 
is observed to correspond with the strong stochasticity 
threshold that exists for non-integrable Hamiltonian sys- 
tems with many degrees of freedom. This transition, sig- 
naled by a crossover in the e-dependence of the largest 
Lyapunov exponent, occurs between dynamical regimes 
which exhibit qualitatively different efficiency of phase 
space mixing. 

The dynamical approach to homopolymer dynamics 
also demonstrates a clear change in the scaling of the ra- 
dius of gyration with N as the self-avoidness parameter 
in the potential was considered in addition to the simple 
nearest neighbor potential that maintains connectivity. 
It was seen that the scaling law depends on the relative 
magnitude of the effective range of the repulsive potential 
and of the equilibrium monomer spacing, which is phys- 
ically expected. Previous methods to determine scaling 
relied on computing the radius of gyration for various 
particle numbers, and thus the dynamical approach pro- 
vides and alternative method by which to study scaling 
relations. 

The dynamical method applied to homopolymers also 
provides a clear signature of the 0-transition between 
filamentary and globular configurations. These results 
are obtained for relatively small N, while the statistical 
mechanical approach relies on the thermodynamic limit, 
N — > oo. The results are physically relevant, especially 
for polymers containing small numbers of monomer units. 

We anticipate that the extension of this paper's tech- 
niques of Hamiltonian dynamics to more complex het- 
eropolymer systems such as minimalist protein models 
(see, e.g., pTl llfij ) will provide useful information to 
complement current statistical mechanical and simula- 
tion studies. The dynamical signature of a transition 
from fast to slow phase space mixing in a homopolymer 
model suggests that a similar dynamical signature indi- 
cating the transition from strong to weak chaos may exist 
in protein-like systems. 
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APPENDIX A: STATISTICAL AVERAGE OF (II 2 ) 
FOR A RWM 



Statistical averages are computed in the canonical en- 
semble, giving the partition function for the system in d 
dimensions as 



. n . N 

Z= / n d P' / II d * exp[-/3ff(q,p)] 



(Al) 



where the Hamiltonian is given by Eq. (yi with potential 
energy ©. The integration over the momenta is straight- 
forward and gives a factor 



(A2) 



The remaining integration is simplified upon a change in 
the variables from the {q^} to the {r^} = {rj,j_i} (the 
Jacobian of this transformation is (N + so that 



Z q = {N + 1) 



2ni 



N 



T(d/2) / 

where T is the Eulcr function and 



(A3) 



Ta = 



r exp 



— ^-(r-r ) -—{r-r 



(A4) 
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The square of the distance from one end of the polymer 
chain to the other can be written 



N 



R 2 = |qw -q | 2 = J2 ri " r J 



(A5) 



Hence for the mean value of the end-to-end distance we 
find 



(R 2 ) = 



N 



IlilxdvieM-PH) 



(A6) 



Not surprisingly, the power law observed in the RWM on 
a lattice is recovered. The coefficient £(/3) represents the 
mean distance between particles in the chain (counter- 
part to the lattice constant in the discrete model). By 
evaluating the integral Eq. <|A4|> . we easily see that l{(3) 
is a monotonically decreasing function of (3, going to ro 
in the limit (3 — * oo. 



APPENDIX B: DEFINITION OF THE LARGEST 
LYAPUNOV EXPONENT 

Given that x l = X 1 (x 1 . . . x n ) is a dynamical system, 
and if we denote by = J^[x(t)]^ k the usual tangent 
dynamics equation, where [J7^] is the Jacobian matrix of 
[X 1 ], then the largest Lyapunov exponent Ai is defined 
by 



t-oo t ||£(0)|| 



(Bl) 



By setting A[z(*), £(t)] = £ T J[x(t)} £/ £ T £ = = 
ln(^ T ^), this can be formally expressed as a time av- 
erage 



Ai 



1 



lim - / dr A[x(t),£(t)} 



(B2) 



APPENDIX C: RIEMANNIAN THEORY OF 
HAMILTONIAN CHAOS 

1. General framework 

We provide a brief summary of the concepts involved 
in a differential geometric method to describe Hamilto- 
nian chaos and to theoretically predict the value of the 
largest Lyapunov exponent. The geometrical formula- 
tion of the dynamics of conservative systems was first 
used by Krylov in his studies of the dynamical founda- 
tions of statistical mechanics 0] and subsequently be- 
came a standard tool in the study of abstract systems 



like Anosov flows in the framework of ergodic theory. In 
recent papers j33, 0| , the geometric approach has been 
successfully extended and applied to explain the origin 
of chaos in models of physical relevance. 

Consider a system with N degrees of freedom defined 
by the Lagrangian C = T — V, in which the kinetic energy 
is quadratic in the velocities, T = ^aijX i x^ . Such dynam- 
ics can be recast in geometrical terms since the natural 
motions are the extrema of either the Hamiltonian ac- 
tion functional, 6>h = / £di, or the Maupertuis action, 
5m = 2 J T At. In fact, the geodesies of a Riemannian 
manifold are themselves the extrema of the arc-length 
functional i = J y/gijdx^dx^. Hence a suitable choice of 
the metric tensor allows identification of the arc length 
with either 5h or 5m, and of the geodesies with the nat- 
ural motions of the dynamical system. 

Starting from 5m we obtain the Jacobi metric on the 
accessible configuration space, {gj)ij = [E— V({x})] a%j. 
The extrema of the Hamiltonian action <Sh can be de- 
scribed as geodesies of a manifold using Eisenhart's 
metric on an enlarged configuration space-time, {t = 
x°, x 1 , . . . , x , x N+1 } 7 where real coordinate x N+1 is re- 
lated to the action |31j . The arc- length is 



ds 2 



-2V(x)(dx° 



aijdx t dx 3 



2dx°dx N+1 . (CI) 



The manifold has a Lorentzian structure and the dy- 
namical trajectories are geodesies satisfying the condi- 
tion ds 2 — Cdt 2 , where C is a positive constant. In 
the geometrical framework, the stability of trajectories 
is mapped to the stability of geodesies, and is thus com- 
pletely determined by the curvature properties of the un- 
derlying manifold. The field J, commonly known as the 
Jacobi field, obeys Jacobi equation 



V 2 J + R(^, J) 7 = 



(C2) 



where V-y is the covariant derivative along the geodesic 
7(s), R is the Riemann curvature tensor, and 7 is the 
velocity field along 7. Thus J is seen to measure the 
deviation between nearby geodesies. 

In the case of isotropic, hyperbolic manifolds, the cur- 
vature term in Eq. I|C2[1 can be rewritten as i?(7, J)j — 
KJ, where the sectional curvature if is a negative con- 
stant, thus giving that equation l|C2l) has exponentially 
growing solutions and that the system is dynamically un- 
stable. This is the origin of chaotic dynamics in Anosov 
flows. The geometric picture for coupled nonlinear os- 
cillators, however, is much different since all curvatures 
(Ricci, scalar, sectional) are mainly (and in some cases 
strictly) positive. Even so, exponentially growing so- 
lutions of the stability equation l|C2() can be obtained 
through parametric resonance even if no negative cur- 
vature is experienced by the geodesies [jnl EH HE EH • 
In the large N limit and under the assumption that the 
manifold is nearly isotropic, this mechanism can be mod- 
eled by replacing Eq. i|C2(l with an effective scalar Jacobi 
equation [33, [34J which reads 



1^ 



k(s) ip = 



(C3) 
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where ip 2 cx | J| 2 , and the "effective curvature" k(s) is a 
stochastic, (5-correlated gaussian process. The mean k 
and variance a K of k(s) are identified respectively with 
the average and the root-mean-square (RMS) fluctua- 
tions of the Ricci curvatures at any given point ItR — 
Kr/N (which is itself the average of the sectional curva- 
ture over the directions of J) along a geodesic, 



N 



((K R -(K R )) 2 ) 
N 



(C4) 



Using Eisenhart's metric K R = AV = d 2 V/dx 2 , 
the exponential growth rate Ai of the envelope of solu- 
tions to Eq. l|C3fl . which provides a natural estimate of 
the Lyapunov exponent, can be computed exactly, with 



Here 2r = (7t^/k^)/(2^k (k + a K ) 
limit (J k /kq <C 1, one finds A cx a 2 . 
found in Refs. [allli- 



+ 7r<7«), and in the 
The details can be 



2. Lyapunov exponent for an ideal homopolymer 

We use Eisenhart's metric to compute the mean and 
RMS fluctuations of the Ricci curvature, 



Kb 



h dxV 

i— ( 



(C6) 



A 2k q 
X 2 ~ 3A 



A = 



2g 2 t- 



64K „ A 9 

- -J- 4air 2 
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(C5) 



which are necessary to determine the largest Lyapunov 
exponent, Ai. In d dimensions, the Hamiltonian de- 
scribed by Eq. Q with potential energy 10 gives 



Kb 



EE 



d 2 V(q) 



j=0 n=l **<n 
N , 

2 E ( a + 36 ( r w-i - r °) 2 + ( d - 1 ) 



a(rjj-i - r ) + b(r j ^ J ^ 1 - r f 



In order to compute the Gibbs average, we modify the The factorization of the partition function along the lines 



canonical partition function Eq. IjAlJ) such that 



of Add. IXI yields 



- N . N 

Z(a) = / Y[d Pl / Y[d qj exp[-/3H(p,q)+aK B ] . 
J i=o •* j=0 

(C7) 

The problem of computing the mean and RMS values of 
the curvature thus reduces to taking derivatives of the 
modified partition function (|C7|I . since 



N 



1 

N 
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In Z(a) 



(C8) 



Q = 



and 



((K R -(K R )) 2 ) 1 



N 



N 



-^-zlnZ(a) 



Q=0 



Z{a) 



(N + l) 



r(d/2) J 



(C9) where 
I 



(CIO) 



id{a) 



+00 



drr d 



exp 



a/3 



(r - r ) 2 



63 



(r — ro) 4 + afc(r) 



(Cll) 



and Thus the quantities to be computed are the one dimen- 

rn ( r _ rn \3 sional integral in Eq. (101 H and its derivative with re- 

jfe(r) = 2ad-2a(d-l) — + 6b(r-r ) 2 + 2b(d-l)± 

r r 

(C12) 
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spect to a, each of which is readily computed numerically. 



